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Abstract 


There  has  been  frequent  controversy  over  the  years  regarding  the  use  of  nu¬ 
merical  rootfinding  for  the  solution  of  queueing  problems.  It  has  been  said 
that  such  problems  quite  often  present  fairly  typical  computational  difficul¬ 
ties.  However,  it  turns  out  that  rootfinding  in  queueing  is  so  well  structured 
that  problems  do  not  occur.  There  are  fundamental  properties  possessed  by 
the  well-known  queueing  models  that  eliminate  classical  rootfinding  prob¬ 
lems.  Most  importantly,  we  show  that  uniqueness  of  roots  is  standard  within 
simply  determined  regions  in  the  complex  plane  and  prove  that  the  G/ Ek/\ 
model  has  unique  roots  easily  found  inside  the  complex  unit  circle.  Extensive 
computational  results  are  given  to  support  our  contentions.  ^ 
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1  INTRODUCTION 

In  all  the  major  textbooks  on  queueing  theory,  we  see  repeated  ref»rence 
the  role  of  rootfinding  in  the  solution  of  many  of  the  primary  models.  A 
large  number  of  these  models  can  be  solved  by  other  numerically-intensive 
procedures  which,  however,  preserve  a  much  closer  link  to  the  underlying 
stochastic  processes.  (See  especially  the  work  of  Neuts,  e.g.  1981.)  The  point 
is  often  made  that  the  use  of  transforms  and/or  generating  functions  leads 
away  from  probabilistic  arguments  to  the  application  of  analytic  function 
theory.  But  much  of  the  ”anti-transform”  position  is  directed  at  the  typical 
need  to  find  roots  to  a  transcendental  or  high-degree  polynomial  equation 
in  order  to  finish  off  a  transform  solution.  It  is  the  common  wisdom  that 
such  rootfinding  is  fraught  with  critical  obstacles,  such  as  difficulties  raised 
by  multiple  roots  or  sharp  slopes  in  the  function,  so  that  some  roots  may 
indeed  not  be  found  by  even  the  most  sophisticated  algorithms.  However, 
it  turns  out  instead  that  rootfinding  in  queueing  is  so  well  structured  that 
these  problems  do  not  occur.  It  is  the  objective  of  this  work  to  show  that  the 
standard  rootfinding  problems  found  in  queueing  analyses  have  fundamental 
properties  that  always  allow  computational  solution. 

It  is  important  to  recognize  that  the  efficient  finding  of  roots  is  essentially 
a  three-phase  process.  First,  we  need  to  determine  their  multiplicity,  that 
is,  whether  or  not  there  are  any  repetitions.  Second,  roughly,  where  are  the 
roots?  (It  is,  for  example,  quite  useful  to  know  how  many  are  inside  or  on 
the  unit  circle  in  the  complex  plane.)  And,  finally,  we  need  to  develop  an 
effective  algorithm  and  employ  it  on  appropriate  hardware  to  obtain  complete 
answers  within  reasonable  amounts  of  time. 

2  SOME  KNOWN  AND  BASIC  RESULTS 

It  has  already  been  shown  in  a  number  of  queueing  models  that  roots  arc 
unique  and  that  they  can  easily  be  obtained  numerically.  (See  Chaudhry 
and  Templeton,  1983,  Chaudhry,  Madill  and  Briere,  1987,  and  Briere  and 
Chaudhry,  1987,  for  example.)  We  present  a  short  discussion  of  some  of  the 
more  familiar  ones  here  (in  no  particular  order)  to  illustrate  what  is  already 
known  and  to  lay  some  groundwork  (with  proofs  streamlined)  for  the  more 
generic  results  which  follow.  Section  3  will  provide  more  complete  details  on 
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root  calculation  for  the  critical  models  under  the  most  extreme  of  conditions. 


2.1  M/M(a'V  1 

Just  about  the  simplest  Markovian  queue  for  which  rootfinding  becomes  ger¬ 
mane  is  the  bulk  variation  of  the  classic  M/M/1  in  which  service  is  in  fixed 
batches  of  size  K  (whether  or  not  the  server  has  to  wait  for  a  full  batch 
of  size  K  is  irrelevant.  The  key  equation  here  appears  in  the  denominator 
of  the  queue-length  generating  function  (or,  equivalently,  as  the  character¬ 
istic  equation  when  the  stationary  system-size  probabilities  are  related  in 
difference-equation  form): 

pzK+1  -  (p  +  l)zK -f  1  =  0  (p  =  A/ p).  (1) 

First,  we  can  show  under  ergodicity  (using  Rouche’s  theorem  -  see  Gross 
and  Harris,  1985,  for  example)  that  this  equation  has  one  root  of  1.  K-l  roots 
inside  or  on  the  unit  circle,  and  finally  one  outside  the  unit  circle.  This  latter 
root  ultimately  determines  the  values  of  the  stationary  probabilities. 

Though  we  do  not  need  to  determine  the  roots  inside  or  on  the  unit 
circle  to  solve  the  queueing  problem,  it  is  useful  to  look  at  them  a  little 
more  closely.  It  turns  out  that  the  roots  are  unique  (as  first  noted  by  Bailey, 
1954),  a  property  quite  typical  of  the  kinds  of  rootfinding  problems  arising  in 
stochastic  models.  This  is  proven  by  contradiction,  and  we  start  by  assuming 
that  z,  ^  0  is  a  repeater.  Since  the  derivative  must  then  vanish  at  z,,  it  follows 
that 


ztK-,[p(A'+  l)zt-(p+l)A']  =  0, 


thus  implying  that 

K  p+1 

z.  = - - — . 

K  +  1  p 

By  assumption,  we  see  that  z,  is  a  repeated  root,  and  it  must  be  inside 
or  on  the  unit  circle.  But 


|  z,  |<  1  iff  p/K  >  1 . 


(2) 
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However,  this  violates  the  condition  for  ergodicity.  Hence  the  roots  are  dis¬ 
tinct.  (Of  course,  we  recognize  that  this  is  the  same  construct  as  that  for  the 
Ek/M/1  queue.) 

In  order  to  find  the  root  outside  the  unit  circle  in  this  case,  you  can 
observe  that  it  must  be  real,  since  it  is  unique.  Moreover,  the  function  is 
0  when  z  =  1  and  the  derivative  of  the  function  is  negative  in  the  interval 
[l,z*],  where  2*  =  A'(p4  1  )/p(K  4  1).  Hence  we  need  examine  only  the  real 
interval  greater  than  2*. 

Newton’s  method  finds  successively  better  approximations  for  the  real 
root  of  F(z)  =  0  via  the  equation  Z(n+i)  =  Z(n)  -  F(z(n))/.F'(z(n)).  A  sufficient 
condition  for  convergence  of  this  method  is  that  the  first  two  derivatives  be 
positive.  This  is  easily  shown  for  our  case  when  z  >  2*: 

F'(2)  =  2K-1[p(A'  +  l)2-(p+l)X] 
has  both  factors  positive;  and 

F"(z)  =  KzK-2[r(K  4  l)z  -  (p  4  1)(K  4  1)] 

is  positive  when  z  >  (p4  1)(A'  -  l)/p(A'4  1).  But  this  ratio  is  smaller  than  z ’ 
and  hence  the  second-order  condition  for  convergence  is  also  satisfied.  Thus 
it  appears  that  Newton’s  method,  perhaps  with  a  starting  value  of  (p4  1  )/p, 
is  the  fastest  and  easiest  way  to  find  the  root  outside  the  unit  circle.  Note 
that  we  wish  to  avoid  2*  as  a  starting  point  because  the  derivative  vanishes 
there.  (Since  the  function  is  a  polynomial,  with  a  second  derivative  easily 
computed,  a  variant  of  Muller’s  method  of  approximation  by  a  quadratic 
might  converge  even  faster.) 

A  second  method  for  finding  the  real  root  outside  the  unit  circle  is  a  fixed- 
point  iteration  or  successive-substitution  approach.  The  original  polynomial 
expression  is  rewritten  as  the  ( K  4  l)st  root  of  [(p  4  1)2^  —  1  \/p  and  then  z 
is  repeatedly  substituted  into  this  expression.  A  sufficient  condition  for  the 
convergence  of  this  method  is  that  the  derivative  of  this  ( K  4  l)st  root  be 
between  -1  and  1,  so  that  the  mapping  truly  contracts  the  domain. 

This  condition  can  be  proven  as  follows: 


1 

(P+  l)z*  -  1 

**  (p  4  l)A'z* 

K  4  1 

-  p 

P 
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K(P  +  l)  K-\(P  ±  l)z* 

(K  +  l)p  l  p 


L^-iC 


We  then  recognize  the  first  factor  as  z*  defined  earlier  and  restrict  our  atten¬ 
tion  only  to  the  interval  where  z  >  z*: 


/'(*)  =  * 


-iKP  +  1  )zK  -  1 


so  that 


/'(*)  <  ^ 


(p  4  l)zK  -  1 


pz 


<  1. 


( p  4  \)zK  -  1 

Thus,  we  see  that  the  condition  is  satisfied  and  this  approach  is  also  guaran¬ 
teed  to  converge  as  long  as  the  starting  value  exceeds  z*. 


2.2  A/W/A//1  or  M/EK/ 1 

Here  the  characteristic  equation  whose  roots  we  need  to  find  is  given  by 


pzK  +  1  -  (p4  l)z  4  1  -  0.  (3) 

(Of  course,  we  recognize  that  the  state  probabilities  can  be  obtained  without 
rootfinding  by  going  directly  to  the  usual  recurrence  relationship  building 
up  from  the  easily  computed  p0,  as  in  Gross  and  Harris,  1985,  Chapter  3  ) 
Equation  3  is  known  to  have  a  root  of  1  and  to  have  the  remaining  K  all 
outside  the  unit  circle.  As  in  Section  2.1,  the  lack  of  multiplicity  can  be 
verified  here. 

First,  note  that  since  any  repeat  root  must  make  the  deriviative  of  the 
left-hand  side  of  (3)  vanish,  it  follows  that  z,  would  be  a  repeater  if 


P  +  1 

P(K  4  1)' 


To  determine  whether  or  not  z,  is  a  multiple  root,  we  need  then  to  see  if  it 
satisfies  (3): 
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0  =  pz*+1  ~  (p  +  l)zt  +  1  =  Z,~~r  -  (p  +  l)zt  +  1 

A  t  i 
or 

?  K  +  1 

2‘~  K{P  +  iy 


Thus  we  need  to  show  that 


K  +  1 


v< 


K(p+1) 


P  +  1 

p(K  +  iy 


or  that 


K  +  1 


Kp. 


But  these  two  terms  cannot  be  equal  since 


a'(p+i)T+1  r 

j  1  -Kp\ 

_  K  +  i  J  [ 

K  +  lj 

(K  +  l)(l  -  Kp) 
K  +  1 


Kp. 


Thus  repetition  is  not  possible  and  the  roots  are  simple. 

The  actual  estimation  of  th”  roots  of  (3)  (or  more  appropriately,  (3)  trans¬ 
formed  by  u  =  l/z  to  get  the  roots  inside  the  unit  circle)  can  be  done  nicely 
using  the  Jenkins-Traub  (1970)  rootfinding  algorithm.  This  is  a  cubically 
convergent  application  of  Newton’s  method  to  a  rational  approximation  of 
the  polynomial.  Zeros  are  calculated  in  roughly  increasing  order  of  moduli, 
which  are  all  less  than  1  on  the  inside  of  the  unit  circle.  This  approach  works 
particularly  well  because  of  the  isolation  of  the  roots  and  the  fact  that  defla¬ 
tion  once  again  leaves  a  polynomial  to  solve.  Purification  would  be  advisable 
for  large  values  of  K  to  make  the  technique  work  more  smoothly. 
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2.3 


M/D/c 

As  noted  in  Gross  and  Harris  (1985),  the  fundamental  root  equation  for 
the  M/D/c  (or  equivalently,  the  Al/D^/1)  queue  is  (with  the  service  rate 
assumed  to  be  1) 


1  -  zcexp[X(l  —  z)]  =  0,  (4) 

or 

zc  =  esrp[-A(l  -  z)].  (5) 

There  are  c  roots  inside  and  on  |  z  |=  1,  one  equal  to  1  and  the  remaining 
c-1  inside  (see  Chaudhry  and  Templeton,  1983).  To  show  that  these  c  roots 
are  again  distinct,  we  create  a  contradiction  beginning  from  the  assumption 
that  z,  ^  0  is  a  repeater  Then  the  derivative  of  the  left-hand  side  of  (4) 
should  vanish  at  z,.  The  derivative  is  easily  seen  to  be 

«_,e*p[A(l  -  z,)])(Az,  -  c), 

which  can  only  be  0  if  z,  =  c/A.  But  this  cannot  happen  since  c/A  is  the 
reciprocal  of  the  system  utilization,  which  must  be  less  than  1  for  ergodicity, 
and  therefore  z,  is  outside  the  unit  circle.  So  the  roots  inside  and  on  the  unit 
circle  are  simple  again. 

To  find  the  roots,  the  following  technique  was  originally  developed  by 
Downton  (1955),  expanded  by  Powell  (1985),  and  later  improved  systemat¬ 
ically  by  Chaudhry  together  with  various  collaborators  in  1986  and  1987. 
Equation  5  is  clearly  equivalent  for  all  n  =  3,2,  ...,c  to 

z  —  exp[-A(l  -  z)/c}exp(2‘nni/c),  (6) 

where  t  is  the  squareroot  of  -1  and  exp(2-nni/c),n  —  1,2 ,  ...c,  are  the  c 
complex  roots  of  unity.  Chaudhry,  Madill  and  Briere  (1987)  have  shown 
that,  for  each  n,  (6)  has  exactly  one  root  inside  the  unit  circle.  Powell  (1985) 
used  Newton’s  method  in  his  work  to  find  that  root;  but  Chaudhry  et  al. 
have  had  more  (and  complete)  success  using  Muller’s  method,  never  having 
encountered  any  numerical  difficulties.  We  present  more  on  a  similar  problem 
later. 

Further  results  are  available  when  there  is  random,  bulk  input,  i.e.,  when 
the  model  is  A/T'  1/D/c.  Equation  5  must  now  incorporate  information  on 
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the  probabilistic  nature  of  the  batch-input  sizes.  This  is  done  through  their 
probability  generating  function,  say  X(z).  Then  the  characteristic  equation 
becomes 


zr  =  exp(-A[l  -  A'(z)]), 

with  c  roots  found  to  be  inside  and  on  the  unit  circle. 

2.4  G/M/l  and  G / EK / 1  or  G^/M/l 

The  waiting-time  distribution  function  for  the  general  single-server  problem 
with  exponential  service  requires  the  lone  real  root  on  (0,1)  for  the  (funda¬ 
mental  branching  process)  equation 

OO 

2  =  YLb'z'  =  '4*l/x0  -  z)i  - /?(*) 

»-o 

where  .4*  is  the  Laplace-Stieltjes  transform  (LST)  of  the  interarrival  times 
and  0  is  the  probability  generating  function  (pgf)  of  the  arrivals  during  ser¬ 
vice  (see  Gross  and  Harris,  1985).  The  LST  is  easily  shown  to  be  monotone 
nondecreasing  and  convex,  and  thus  the  root  is  readily  obtainable.  For  exam¬ 
ple,  the  well-known  Newton-Raphson  approximation  method  is  guaranteed 
to  converge  because  of  the  convexity. 

The  problem  becomes  more  interesting  when  Erlang(K)  service  times  are 
used  instead.  Here  the  roots  need  to  be  generally  located  and  then  found  for 

OO 

zK  =  ]T6,z  -  A*\n{\  -  z)j  =  0(z)  (7) 

•=o 

or,  using  z  —  rexp(-id)  and  principal  values  for  all  log  calculations, 

Kln(r)  -f  iK6  =  ln{A'{p{  1  -  i  jzp(i#))]}  +  2tttv  (8) 

for  n  —  (see  Chaudhry,  Madill  and  Briere,  1987).  There  is  clearly 

a  root  at  unity,  and  by  Rouche’s  theorem,  we  can  once  again  show  that  there 
are  K  others  inside  the  unit  circle  |  z  |=  1.  For  each  n  now,  Chaudhry, 
Jain  and  Templeton  (1987)  note  th "* t  there  is  a  unique  root  with  absolute 
value  less  than  1,  using  the  complex  version  of  the  monotonicity  and  convexity 
argument  employed  for  the  G/M/l. 
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Note  that  the  characteristic  equation  for  the  model  G^x^/M/ 1  bears  a 
great  similarity  to  (7).  We  know  that  zK  is  the  probability  generating  func¬ 
tion  of  the  batch-input-size  distribution  since  the  only  possible  size  is  iden¬ 
tically  K .  Rewrite  (7)  then  as 


l  =  0(z)(l/z)K. 

Then  replace  the  term  (1  / z)K  by  the  input-size  pgf  X(z)  evaluated  at  1/z, 
thus  giving  the  c.e. 


1  =(3(z)X(l/z) 

The  roots  of  (8)  are  found  by  separately  solving  its  real  and  imaginary 
portions.  We  know  that  there  is  a  unique  answer  when  (8)  is  evaluated  for 
individual  values  of  n.  But  prior  to  this  work,  it  was  assumed  that  roots 
obtained  in  this  manner  could  be  the  same.  We  have,  however,  now  been 
able  to  show  that  repetition  is  not  possible!  We  provide  our  proof  of  this 
assertion  in  Section  3.2. 

Partial  results  on  the  uniqueness  of  all  the  roots  for  this  model  type  are 
available.  It  is,  in  fact,  known  that  the  roots  are  totally  unique  when  the 
interarrival  times  are  Erlang(J)  distributed  (see  Chaudhry  and  Templeton, 
1983).  In  this  case,  Equation  7  becomes 

zK  =  \Jp/{Jp+K(\-z)}J.  (9) 

Again,  we  need  to  have  the  derivative  vanish  at  any  repeated  root,  say  z,. 
This  is  equivalent  to  requiring  that 

*.K_l  =  (if) 

When  (9)  is  divided  by  (10),  we  see  that 

z,  =  (J p  -I-  h  )/(J  -I-  h). 

When  this  value  is  substituted  back  into  Equation  9,  we  find  that  the  condi¬ 
tion  for  repetition  is  equivalent  to  requiring  that 

Pc  ~  pc  +  (1  -  c).  (11) 
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The  right-hand  side  of  (11)  is  a  straight  line  in  p,  with  y-intercept  of  1-c 
and  slope  of  c,  while  the  left-hand  side  is  a  simple  monomial  with  integral 
power  c.  The  two  functions  intersect  only  at  p  —  1,  which  would  violate  the 
condition  for  ergodicity,  so  that  the  assumption  of  repetition  must  be  false. 

But  not  much  more  has  been  generally  known  about  the  effect  of  the 
form  of  the  interarrival  distribution?  Chaudhry,  Jain  and  Templeton  (1987) 
tried  the  method  of  Equation  8  on  a  large  number  of  cases  and  did  not  en¬ 
counter  any  problems.  They  used  the  QPACK  software  package  developed 
by  Chaudhry  and  his  collaborators  (Chaudhry,  1988).  However,  such  com¬ 
putational  testing  clearly  cannot  be  exhaustive.  We  address  this  question  in 
a  more  complete  fashion  in  Section  3. 

The  waiting-time  distribution  function  for  the  G/Er/I  system  (see,  for 
example,  Chaudhry  and  Templeton,  1983)  is 

K 

VV,(<)  =  1  -  C]T  a,rtea;p[-/z(l  -  r,)t]/(l  -  r.) 

.=  i 

where  C  is  the  arrival-point  probability  that  there  are  no  phases  present,  and 
the  {r,}  are  the  K  complex  roots  inside  the  unit  circle.  The  mixing  constants 
{a,,t  =  1,2,..., r}  are  found  by  the  formula 

a,  =  rW(r‘  -rr)  (r**rr)- 
2.5  M/E(P/ 1,  M/G(y)/\  and  EK/G/l 

For  the  more  general  batch-service  system  with  Poisson  input,  Erlang(K) 
distributed  service  times  and  random  batch  sizes,  the  characteristic  equation 
whose  roots  we  require  is 


(p  +  X  -  Xz) 


^r(i/2)  -1=0, 


(12) 


where  Y(z)  is  the  probability  generating  function  for  the  random  service 
batch  (defined  to  have  finite  support).  If  we  write  Y(z)  as  the  polynomial 


Y(z)  =  y\z  y2z 2  +  y3z3  +  ...  +  yszfc, 
then  Equation  12  can  be  rewritten  as 


9 


* 


Zb  =  [l+  -(1  -  ^)j-,f(y^^b-1  +  y2zb~2  +  ...  +  yb). 

A* 

Clearly,  1  is  a  root;  in  addition,  K  roots  are  outside  the  unit  circle,  while 
b-1  are  inside  and  on.  All  of  the  roots  inside  are  again  distinct  and  easy  to 
obtain  (see,  for  example,  Chaudhry  and  Templeton,  1983).  But,  as  noted 
in  our  A1  /I  discussion,  it  is  the  roots  outside  which  are  critical.  The 

model  A//(7(y  tyl  is  the  generalization  of  the  E k  service  model  and  has  any 
number  of  roots  outside  the  unit  circle  (e.g.,  for  G  =  M,  there  is  1). 

Briere  and  Chaudhry  (1989)  have  analyzed  in  detail  relatives  of  this  model 
where  service  is  instead  either  hyperexponential,  uniform  or  constant.  Re¬ 
sults  are  once  more  quite  favorable. 

To  show  the  further  equivalence  of  this  model  type  to  the  Ek/G/\,  let 
1'  =  K  in  the  M/G ^  V 1  and  then  convert  the  constant  batch  size  and  Pois¬ 
son  input  over  to  an  Erlang(K)  input.  The  resultant  characteristic  equation 
is  thus 


=  B’|A(1  -  *)],  (13) 

w.iere  B*  is  the  Laplace-Stieltjes  transform  of  the  service-time  distribution. 
Note  the  clear  analogy  to  the  c.e.  for  G/Ek/  1  given  by  Equation  7.  Equa¬ 
tion  13  (with  K  =  b)  is  also  the  characteristic  equation  of  the  model  M /G0|t>/1, 
where  the  server  takes  batches  of  size  b  if  available  and  otherwise  waits  until 
at  least  a  customers  are  waiting. 

It  is  well  known  that  (13)  has  K  roots  inside  and  on  the  unit  circle, 
including  the  root  z-1  (see,  for  example,  Abolnikov  and  Dukhovny,  1987). 
We  show  here  that  z  =  l  is  simple  using  the  usual  derivative  test.  To  do  so, 
we  evaluate 


KzK~]  =  - -  \z)ldz 

at  z=l  and  find  that 

K  =  —  A(-l /n)  =  Kp, 

or  p  —  1,  which  is  a  contradiction.  Hence  z=l  cannot  be  a  double  root. 
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When  the  service  times  are  deterministic,  we  can  further  show  that  all 
K-l  roots  inside  or  on  the  unit  circle  are,  in  fact,  strictly  within.  This  follows 
when  we  rewrite  (13)  as 


ZK  _  eKp{z- 1) 

and  assume  that  z,  has  absolute  value  of  1  but  is  not  precisely  equal  to  1. 
Then  we  see  that 


1  =  eK>(*.-i), 

which  implies  that  Re(zt)  =  0  and  thus  that  z,  =  1.  But  this  is  contrary  to 
our  earlier  verification  that  the  root  z  =  1  is  simple. 


3  FURTHER  RESULTS  AND  COMPUTA¬ 
TIONAL  EXPERIENCES 

In  this  section,  we  provide  what  we  feel  is  (fairly)  complete  evidence  on  the 
ultimate  efficacy  of  rootfinding  in  queueing.  For  each  of  the  models  of  Section 
2,  we  shall  expand  on  the  theory  of  root  location  and  repetition,  and  then 
computationally  push  each  problem  to  its  extremes  to  reinforce  the  notion 
that  roots  can  always  be  found  for  queueing  models. 

3.1  A/W/A//1 

First,  we  report  on  further  computational  experience  with  the  Markovian 
batch-arrival  model.  (We  do  not  do  likewise  for  the  Markovian  batch-service 
queue  since  it  requires  the  location  of  only  one  root  and  that  is  not  difficult.) 
We  are  thus  working  with  the  polynomial  defined  in  Equation  3,  and  have 
used  the  software  developed  by  Chaudhry  and  Hasham  (1987).  The  key 
parameters  are  thus  the  bulk  size  K  and  the  traffic  intensity.  We  have  selected 
the  extreme  values  of  the  traffic  intensity  for  our  experiments  to  be  .05,  .1, 
.9  end  .95.  We  have  also  run  at  the  value  .5  as  a  calibration.  The  values 
of  K  chosen  were  10,  25,  50  and  100.  The  results  are  presented  below  in 
Table  1  for  runs  performed  on  an  80286-based,  12  mHz  AT  clone  with  a 
math  coprocessor.  Here,  as  for  subsequent  examples,  root  values  have  been 
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verified  to  be  correct  using  the  IMSL  package  ZANLY  for  analytic  function 
rootfinding. 


Table  1 

/  M  / 1  Computations 


Intensity 

Batch  Size 

Run  Time  (mm:sec) 
for  Roots 

.50 

10 

0:16 

25 

0:21 

50 

0:28 

100 

0:40 

.05 

10 

2:03 

25 

2:15 

50 

2:25 

100 

2:41 

.10 

10 

1:04 

25 

1:12 

50 

1:19 

100 

1:34 

.90 

10 

0:12 

25 

0:15 

50 

0:22 

100 

0:35 

.95 

10 

0:11 

25 

0:15 

50 

0:22 

100 

0:34 

500 

15:41 

3.2  G/EkI  1 

We  know  that  this  model  has  K  roots  inside  the  unit  circle,  and  as  promised, 
we  now  show  that  these  values  are  indeed  unique. 

Theorem  1  The  roots  of  the  characteristic  equation  of  the  G/ Ek/1  model 
(or,  equivalently,  the  Gi'K'>  /  M  /  \ )  are  unique,  with  one  real  for  K  odd  and  two 
real  for  K  even. 
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Proof:  Use  (7)  in  the  form  zK  =  /3(z).  Then  by  a  geometric  argument 
essentially  the  same  as  that  for  the  G/M/l  used  in  Figure  5.1  of  Gross  and 
Harris  (1985),  it  follows  that  there  exists  a  unique  real  root  in  (0,1)  for  all 
K  when 


/?'(!)>  \d(zK)/dz\I=i. 


But  this  is  equivalent  to 


p  ,,  KX 

-  >  h  or  -  <  1, 

A  p 


which  is  true  from  ergodicity.  [For  K  even,  it  is  easy  to  see  that  there  is  an 
additional  real  root  in  (-1,0).] 

Next,  remember  from  Section  2.4  that  Equation  8  has  a  unique  root  inside 
the  unit  circle  for  each  n  —  1,...,A';  call  it  ( rn,9n ).  But  it  is  also  true  that 
(8)  has  a  unique  (possibly  non-integer)  value,  n„  for  each  pair  (r,,0,).  Thus 
if  we  assume  for  t  ^  j  that  (r, ,#,)  =  (r},8}),  it  follows  that  n,  =  n}.  But 
this  contradicts  the  uniqueness  of  (r,0)  for  each  n.  Therefore  all  K  pairs  of 
roots  (rn,0„)  must  be  different. 

Because  of  this  uniqueness,  we  see  that  the  waiting-time  distribution  func¬ 
tion  of  Section  2.4  for  the  line  delays  is  a  generalized  hyperexponential.  The 
mixing  constants  {a,}  of  the  CDF  formulation  are  guaranteed  to  be  distinct 
and  easily  computed. 

To  show  the  ease  with  which  G/ EK/\  roots  may  be  found,  we  consider  ex¬ 
amples  using  the  most  general  distribution  classes  found  in  queueing,  namely, 
the  phase  types  of  Neuts  (see  Neuts,  1981),  the  generalized  hyperexponentials 
of  Harris  (see  Botta,  Harris  and  Marchal,  1987)  and  the  Coxian  distributions 
(Cox,  1955).  If  the  results  for  the  most  difficult  problems  using  these  kinds 
of  distributions  (with  each  class  known  to  be  dense  in  the  set  of  all  CDFs) 
are  favorable,  then  we  become  much  more  comfortable  in  rootfinding  efforts 
using  any  other  interarrival  distributions.  Two  examples  of  each  class  have 
been  studied,  and  the  results  are  presented  in  the  following.  The  first  PH 
example  is  taken  from  Botta,  Harris  and  Marchal  and  has  CDF 


F(t)  =  1  -  1  .C93e~4  8464  +  0.343e~4  195(  -  0.05e"O  959(. 
Equation  7  now  becomes 
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zK  =  6.266/[4.846+/i(l-*)]-1.439/[4.195+/i(l-z)]  +  .04795/[.959+/i(l-*)]. 

We  have  chosen  three  medium  to  large  values  of  K,  namely,  10,  15  and  30, 
remembering  that  it  is  K  that  determines  the  number  of  roots  we  need  to 
find.  Each  K  has  been  paired  up  with  a  high  and  low  traffic  intensity.  Recall 
that  we  are  to  determine  K  roots  inside  the  unit  circle.  We  have  used  the 
Chaudhry  QPACK  software,  with  all  of  the  roots  found  quickly  in  each  case. 
Our  experiences  are  given  in  Table  2,  and  a  plot  is  presented  as  Figure  1  of 
the  actual  location  of  the  roots  over  the  unit  circle  for  one  of  the  K  =  15 


examples  ( p 

Batch  Size 

—  .1). 

M 

Table  2 

PH/Ex/ 1  Computations 
Intensity 

Run  Time  (min:sec) 

10 

46.0 

.91522 

for  Roots 

0:21 

10 

421.0 

.10000 

0:21 

15 

70.0 

.90214 

0:30 

15 

631.5 

.10000 

0:28 

30 

140.0 

.90214 

0:49 

30 

1263.0 

.10000 

0:49 
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Figure  1:  Root  Location,  I 
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RADIUS  OF  CIRCLE  :  U80000E+01 


15 

ROOTS 


The  first  GH  illustration  is  an  example  offered  by  Botta,  Harris  and  Mar- 
chal,  which  turns  out  to  be  also  PH,  though  not  a  mixed  generalized  Erlang. 
This  is 


F(t)  =  1  -  6e~4t  +  13e~3'  -  8e~2t, 

with  E\T]  =  7/6.  The  resultant  version  of  (7)  which  is  found  when  the 
appropriate  transforms  are  taken  is 

zK  =  24 / [4  +  /x(l  -  *)]  ~  39/(3  +  n{  1  -  z) ]  +  16/(2  +  /t(l  -  *)]. 

The  problem  is  again  most  challenging  when  K  is  at  least  fairly  large.  So 
once  more  we  have  set  K  in  separate  runs  to  be  10,  15  and  30  and  used  both 
a  high  and  low  traffic  intensity  for  each  value  of  K.  We  have  again  used  the 
Chaudhry  QPACA  software  and  all  roots  were  found  quickly  in  each  case, 
with  the  results  displayed  in  Table  3.  For  illustration,  a  plot  of  the  actual 
location  of  the  roots  over  the  unit  circle  for  the  final  of  these  six  examples 
(with  K  =  30  and  p  =  .10006)  has  been  included  as  Figure  2. 

Table  3 

GH/Ek/1  Computations 


Batch  Size 

Intensity 

Run  Time  (min.sec) 
for  Roots 

10 

9.00 

.95238 

0:21 

10 

85.00 

.10084 

0:21 

15 

14.01 

.91771 

0:52 

15 

128.00 

.10045 

0:29 

30 

29.00 

.88670 

1:30 

30 

257.00 

.10006 

0:48 
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Figure  2  Rr>ol  Location, 


RADIUS  OF  CIRCLE  :  0.100800E+01 
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We  also  experimented  with  a  second  GH  distribution  taken  from  Botta, 
Harris  and  Marchal: 

F(t)  =  1  -  4e-‘  +  6e~2t  -  3e-3t. 

This  particular  one  is  of  potential  special  interest  since  it  turns  out  that  it  is 
not  also  a  phase-type  CDF.  The  mean  interarrival  time  1/A  is  found  to  be  2. 
We  have  again  chosen  a  range  of  medium  to  large  values  of  K,  namely,  10, 
15  and  30,  and  have  used  a  pair  of  high  and  low  traffic  intensities  for  each. 
These  results  are  displayed  in  Table  4,  with  a  ploi  of  the  X  =  30,  p  =  .1  case 
presented  as  Figure  3. 


Table  4 

Second  GH/Ek/1  Example 


Batch  Size 

T 

Intensity 

Run  Time  (min 
for  Roots 

10 

5.60 

.89286 

0:21 

10 

50.00 

.10000 

0:21 

15 

8.15 

.92025 

0:53 

15 

75.00 

.10000 

0:28 

30 

16.66 

.90004 

1:31 

30 

150.00 

.10000 

0:49 
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Figure  3  Root  Location,  III 


RADIUS  OF  CIRCLE  :  0.100000E+01 
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Next,  we  went  back  to  a  phase-type  input  process  and  used  the  harmonic 
distribution 

F(t)  =  1  +  e-2  8846j.3868sm(.5897<)  +  .1729coa(.5897t)) 

-1.1729e"  23071 

as  an  illustration  Results  were  again  excellent,  with  all  roots  found  quickly 
and  efficiently  -  see  Table  5 


Table  5 

Second  Phase-Type  Example 


Batch  Size 

P 

Intensity 

Run  Time  (mtn.sec) 
for  Roots 

10 

20.000 

.10000 

0:25 

10 

2  200 

.90909 

0:26 

15 

30.000 

.10000 

0.35 

15 

3.330 

.90090 

0:38 

30 

60.000 

.10000 

1:03 

30 

6.667 

.89996 

1:06 

The  final  two  examples  of  this  section  are  Coxian  distributions,  with  the 
requisite  rational  LSTs.  Each  of  these  is  also  neither  a  GH  or  PH  distribution. 
The  presence  of  an  harmonic  term  guarantees  that  the  density  cannot  be  GH, 
while  the  fact  that  the  df  hits  the  t-axis  at  at  least  one  point  guarantees  that 
we  do  not  have  a  phase-type  problem  as  well.  The  first  case  here  is  the 
density 


f{t)  =  lOe  21  / j  1  cos(t)}. 


All  of  the  roots  were  once  again  found  quickly  and  efficiently  -  see  Table  6. 
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Table  6 

Second  PH /Ex/ 1  Example 

Batch  Size 

Intensity 

Run  Time  (min.sec) 

for  Roots 

10 

50 

.28571 

0:24 

10 

18 

.79365 

0:24 

15 

189 

.11337 

0:32 

15 

23 

.93168 

0:59 

30 

371 

.11552 

1:35 

«,  • 

48 

.89286 

1:23 

The  conch 

jding 

example  is  a  Coxian 

distribution  function  with  an  atom 

at  the  origin: 

Fit)  -  1  -  |(e 

'~l  +  e  21  )■ 

The  results  follow. 

Table  7 

Second  Coxian  Example 

Batch  Size 

P 

Intensity 

Run  Time  fm:n:sec) 

for  Roots 

10 

200 

.10000 

0:19 

10 

22 

.90909 

0:20 

15 

300 

.10000 

0:26 

15 

32 

.93750 

0:27 

30 

600 

.10000 

0:46 

30 

62 

.96774 

0:45 

As  a  last  test 

case 

for  this  example,  we 

used  a  batch  size  of  100  and  service 

rate  of  220  (giving  a  traffic  intensity  of  .90909).  The  run  took  2  minutes  and 
11  seconds,  and  the  roots  are  displayed  in  Figure  4. 


I 

i 

Figure  4.  Root  Loration,  IV 

RADIUS  OF  CIRCLE  :  0.100000E+01 


3.3  G/GEk/1 

An  obvious  extension  of  the  G/£#/ 1  model  is  the  broader  class  of  models 
where  service  times  are  generalized  Erlang(K),  that  is,  they  are  found  as 
the  convolution  of  K  independent  but  not  identically  distributed  exponential 
random  variables.  For  this  very  large  and  dense  class,  we  can  show  as  in  the 
prior  model  that  there  is  exactly  one  real  root  of  the  characteristic  equation 
for  K  odd  and  two  real  roots  for  K  even.  Unfortunately,  uniqueness  may  not 
obtain  here,  but  root  location  is  not  difficult  anyway.  We  thus  present  the 
following  result  on  real  roots  as  Theorem  2. 

Theorem  2  Tht  characteristic  equation  of  the  G/GEr/I  model  has  a  unique 
positive  real  root  when  K  is  odd  and  two  unique  real  roots  when  K  is  even. 

Proof  :  The  c.e.  for  this  model  is  more  commonly  known  from  the  G/G/ 1 
formulation  as 

A'(~s)B‘(s)=  1,  (14) 

where  A *  and  B'  are  the  Laplace-Stieltjes  transforms  of  the  interarrival  and 
service-time  distributions,  respectively.  From  Gross  and  Harris  (1985),  for 
example,  we  know  that  (14)  has  K  roots  with  negative  real  parts.  Rewrite 
(14)  as 


A*(s) 


1 

b^7) 


(15) 


where  the  generalized  Erlang’s  phase  rates  have  been  placed  in  ascending 
order  as  HuPz >  -.Mk  Then  let  s  =  p/<-(l  —  z)  and  substitute  into  (15).  It 
follows  that 


/,  m  n  Mi  Mk  4"  Pkz 
A  (/iK(l  -  *)  =  II - • 

,=i  M. 

But  we  have  thus  effectively  reduced  the  problem  to  one  which  is  nearly  iden¬ 
tical  to  that  of  the  Gj Ex /l  model  of  Equation  7,  namely,  an  LST  evaluated 
at  jr(l  —  z )  equal  to  a  polynomial  of  degree  K.  The  geometric  argument  of 
Theorem  1  now  goes  through,  with  appropriate  real  solutions  as  long  as  the 
expected  number  of  arrivals  per  service  time  is  less  than  1.  □□  . 
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3.4  Ek/G/1  and  M/Ga’b/1 

Remember  from  (13)  that  these  models  have  characteristic  equation 


zh  =  zK  =  B\X  -  X z), 

where  B*  is  the  Laplace-Stieltjes  transform  of  the  service-time  distribution. 

We  next  offer  a  numerical  example  for  this  model  type.  Consider  an 
Ek/G/1  problem  in  which  service  is  generalized  Erlang.  Let  the  GE  distri¬ 
bution  have  phases  with  mean  rates  jij  =  12,  /r2  =  6  and  /x3  =  4.  Then  the 
equation  of  interest  becomes 

K  i  11 

2  ~  1  +  12A(1  -  z)  1  +6A(1  -  z)  1  +  4A(1  -  *)' 

W'e  have  varied  K  in  our  usual  way,  combined  with  a  variety  of  input  rates 
and  thus  traffic  intensities.  Roots  have  once  more  been  found  quite  easily. 
The  results  are  offered  in  Table  8. 


Batch  Size 

X 

Table  8 

Ek/G/\  Example 
Intensity 

Run  Time  ( min.sec ) 

10 

2 

.1 

for  Roots 

0:23 

10 

18 

.9 

0:22 

15 

3 

.1 

0:31 

15 

27 

.9 

0  31 

30 

6 

.1 

0:51 

30 

54 

.9 

0:50 

4  CLOSING  REMARKS 

We  have  tried  to  show  that  rootfinding  in  classical  queueing  models  is  not  a 
difficult  matter.  Our  view  of  this  is  based  on  the  facts  that  the  critical  roots 
typically  appear  singly  in  easily  located  regions  and  that  our  empirical  nu¬ 
merical  experience  has  always  been  most  favorable.  While  it  is  true  that  we 
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have  not  examined  multi-server  queues  in  great  detail,  rootfinding  in  those 
cases  is  generally  quite  similar  to  the  single- channel  models  with  similar  in¬ 
put  and  service-time-distribution  combinations.  There  are  clearly  problems 
whose  basic  characteristics  are  so  involved  that  rootfinding  is  almost  guar¬ 
anteed  to  be  perverse.  But  we  believe  that  the  models  we  have  explored  are 
very  comprehensive  in  their  application  and  that  lessons  from  their  solution 
are  very  useful  in  the  effective  solution  of  the  more  complicated  models. 
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